Nonuniversality of the dispersion interaction: analytic benchmarks for van der Waals 

energy functionals 
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We highlight the non-universality of the asymptotic behavior of dispersion forces, such that a 
sum of inverse sixth power contributions is often inadequate. We analytically evaluate the cross- 
correlation energy E^ r between two 7r-conjugated layers separated by a large distance D, within the 
electromagnetically non-retarded Random Phase Approximation, via a tight-binding model. For two 
perfect semimetallic graphene sheets at T=0K we find E^ r oc D~ 3 , in contrast to the "insulating" 
D~ 4 dependence predicted by currently accepted approximations. We also treat the case where one 
graphene layer is replaced by a thin metal, a model relevant to the exfoliation of graphite. Our 
general considerations also apply to nanotubes, nanowires and layered metals. 
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Dispersion (van der Waals, vdW) interactions 0,0 are 
especially significant in soft matter, though they coexist 
there with stronger covalent and metallic forces. It is 
known 

H H H SB that all the relevant physics ( in the 
electromagnetically non-retarded regime) is contained in 
microscopic generalizations of the LifshitzQ approach: 
these use the fluctuation-dissipation theorem (FDT) to 
generate the electronic correlation energy, starting from 
a sufficiently accurate nonlocal electronic density-density 
response function x- They are "universal" in that they 
capture the vdW interaction in systems of any size, shape 
or composition, and "seamless" in that they work equally 
well at all separations of the sub-systems. Here we point 
out that many currently popular approaches do not share 
this universality: for example, the usual sum of asymp- 
totic atom-atom contributions of form —CqR^ b (see e.g. 
0) is inappropriate in many anisotropic systems of great 
current interest, even in well-separated regimes. 

To implement the FDT aproach, \ can be usefully 
modelled, for example, by time-dependent density func- 
tional theory^ with its exchange-correlation kernel f xc . 
The simplest case, f xc — 0, gives the Random Phase Ap- 
proximation (RPA) correlation energy. While theories of 
this class have been applied successfully without further 
approximation to vdW interactions in simple systems 
such as small dimers|9l Ho| or layered jellium [jllllLll^|. 
the numerics become formidable for soft systems of re- 
alistic complexity. There has been progress recently in 
creating more tractable theories 0, 0, EI 0, El • Based 
on these approaches, functionals have been proposed [lfij 
that are solely or principally density-based, that give sen- 
sible results for the energetics of compact systems and 
that have the typical asymptotic l/R e form 
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for separated electron densities n A = nA(r A ),n B = 
n B (r B ). These seamless theories are promising for 
weakly bonded finite molecules and small clusters, giving 
a natural "saturation" of the R~ 6 van der Waals energy 
contributions at shorter distances, a long-anticipated 
goal^^- Perturbation theory also supports QJ. 

It has therefore been widely supposed that JU is the 
correct asymptotic form, and that only detailed appli- 
cations remain to be analyzed. We show here that the 
situation is not so simple, at least for soft layered and 
striated structures. We motivate our general point by 
new results for the cohesion energetics of graphitic sys- 
tems fill, a controversial El E| and technologically 
relevant |2ll E3 , topic. Graphene sheets are highly po- 
larizable, so vdW forces should be involved |24|. Never- 
theless, the electronic local density approximation (LDA) 
gives good equilibrium layer spacings 25] and breathing 
phonon frequencies 26] despite the complete lack of dis- 
tant vdW forces in the LDA0,0]. However, for the depth 
— min/j E(D) of the binding energy curve of graphene 
sheets distant D, the LDA gives a value (« 20mH / atom) 
that is little more than half of the most accepted exper- 
imental values (35 mH/ atom or more [Hill). Recent 
calculations [12j found the LDA similarly underestimates 
the binding of a pair of thin parallel high-density jel- 
lium electron gases, another strongly anisotropic layered 
system that might (e.g.) represent alkali-intercalated 
graphite. For other layered jellium systems the LDA 
gave varying degrees of unsatisfactory behavior near 
equilibrium 12], related to the absence of vdW interac- 
tions in LDA. Nevertheless the LDA in fact performs 
much better than the common gradient functionals for 
these systems 0, H^J- Current opinion seems to 
be that LDA should be used near the equilibrium lattice 
configuration, and extended by fitting to a sum including 
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atom- atom terms —CeR A % in order to treat "stretched" 
soft matter. 

A known exception is the interaction energy E of two 
thin metals at large separation D, for which CijR~j 6 
predicts E oc —C^D~ 4 but the correct dependence (e.g. 
from RPA energetics |E IH, ) is -G 5/2 £>- 5 / 2 . We 
will show that this is not an isolated case. Qualitatively 
speaking, part of the reason for this difference is that 
conducting electrons can move long distances and are 
not confined to the vicinity of atoms, in contrast to the 
models that lead to R~ 6 force laws. The low dimension- 
ality is also important, leading to weak metallic screening 
and hence to gapless two-dimensional (2D) plasmons at 
low wavenumber. By contrast, thick metallic slabs have 
gapped surface plasmons and give a conventional attrac- 
tion law, E » —C/D 2 whose form can be obtained by 
summing R~ 6 atom-atom contributions. 

What does this type of physics imply for the energet- 
ics of layered planar graphene-based systems? Firstly, 
an isolated graphene sheet is not a metal but a zero- 
gap insulator: the two-dimensional Bloch electron dis- 
persion relation has a "conical" shape near the Fermi 
level, with a zero energy gap at the Fermi level but a 
zero density of states at the gap[2l| . When the graphene 
layers are brought to their equilibrium separation to cre- 
ate graphite, the graphene bands overlap slightly at the 
Fermi level to create small pockets of electrons and holes. 
Thus the layers in equilibrium graphite are metallic, but 
this only involves a small fraction of the electrons. When 
the layers are separated (at zero electronic temperature) 
this overtly metallic character is lost, so that one cannot 
argue for a —C^jiE^I 2 energetics at large layer sepa- 
ration D and T = OK. We show below, however, that 
the attraction energy between well-separated graphene 
planes at T — OK is of form — C^D^ 3 , closer to metallic 
D~ 5/>2 behavior than to insulating D~ 4 behavior. This is 
a principal quantitative result of the present work, and it 
poses a significant asymptotic constraint on electromag- 
netically non-retarded van der Waals energy functionals. 
Fitted interactions including C^i?" 6 terms 18;}, and 



seamless energy functionals proposed recently [71 llq. for 
example, produce the insulating D~ A behavior, and so do 
not conform to the new constraint. Neither do ground- 
state LDA calculations. We now give the details of our 
argument leading to this new D~ 3 form of interaction. 

All the new physics here comes from electrons close 
to the Fermi level: we can ignore the response of the 
tightly-bound covalent sp 2 electrons, whose finite energy 
gap ensures that they produce a conventional vdW at- 
traction of insulator type (energy cx — -D~ 4 ), negligible 
at large separations compared with the D~ 3 vdW attrac- 
tion that we shall find between the tt z electrons of interest 
here. In the tight binding (TB) model, the valence and 
conduction n z Bloch orbitals and their energies e^ v ' c ' (k) 
are constructed via bonding and antibonding combina- 



tions of a single tt z (Wannier) orbital w(r) = w(x + zk). 
From Rayleigh-Schrodinger perturbation theory the zero- 
temperature density-density response xks of indepen- 
dent 7r z electrons moving in the groundstate Kohn-Sham 
potential of a graphene layer is then of the form 

XKs(q,G,G',z,z',co = iu) = SS'*x (q,iu). (2) 

Here (flics in the 2D Brillouin zone, G. G" arc 2D recipro- 
cal lattice vectors, z and z' are positions measured in the 
direction perpendicular to the planes. S = S(q + G,z), 
S' = S{q + G',z') with S(k,z) = j d 2 xe lL3 w*{x + 
zi)w{x + zi), and the effective 2D response in is an 
integral over the 2D Brillouin zone (BZ), 



Xo(q,iu) = - 



\L(p,p+g)\ 2 2 As d 2 p 
p ~ eBZ (2tt) 2 Ae 2 + (hu) 2 - 



(3) 



Here Ae = e^ c '(p) — e^ v '(p + q), and L is an over- 
lap integral evaluated below. L approaches zero when 
\q\ « \p\ , \p + q\ by valence-conduction orbital orthog- 
onality, so that, for q << G, the Brillouin zone inte- 
gral © is dominated by small p values. (Here p is 
measured from a K point in the Brillouin zone where 
the bands cross). For the small q values required be- 
low, we can therefore use the linear band dispersion 
e( v > c )(p) = =fhv \p\. (We take the graphene velocity to 
be v = 5.7 x 10 5 m/s). 



the overlap element is \L(p, p nl " 



211 In the same approximation 
4sin 2 (^) where 9 is 
the angle between p and p 1 . Eq. 10 can be evaluated 
analytically (see e.g. |13|) by extending the p integration 
to all space, finite-BZ corrections vanishing as q — > 0: 



Xo(<7,«u) 



(4) 



The divergence of the polarizability q~ 2 Xo at low q and 
u, and the parallel dependence on q and u, are crucial. 

For two non-overlapping graphene sheets centered at 
z = and z = D, the "bare" (Kohn-Sham) polarizability 
from (J2J is then 

Xo(q, G,G', z, z', iu) = xo(q, iu) {SqSq + S D S^) (5) 

where S D = S(q + G,z- D) and S' D = S(q + G',z' - 
D). The interacting RPA density-density response func- 
tion xx(q>G,G',z,z',uj) satisfies the following screen- 
ing equation, where stars represent both convolution in 
(z, z') space and matrix multiplication in (G, G') space: 



X\ =Xo + Xo* W coul * X \- 



(6) 



Here V coul = 27re 2 exp(-|G + q\\z - z'\)8 QS , /\G + q\. 
The factorizable (z, z') and (G, G') dependence in (J5J 
arises because only a single ir z orbital is involved in 
the TB description of the 7r-conjugated Bloch states. 
This leads to an analytic solution of ©, where x\ 
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has intra-layer terms xaii, XA22 similar to J5J, and 
also inter-layer terms XA12, X\2i- Xaii = XA22 = 
(1 - Xo^ii)Xo/e, XA12 = XA21 = V 12 x%/e, where 

s(q, iu) = (1 - xo{Vu + V 12 ))(l - MV11 ~ V 12 )), (7) 

ViM = E / dz"dz'"S* (q + G", z" - Zi) 

G " 

xV coul {q + G",z" - z"')S{q + G",z"' - Zj), (8) 

and with Z\ = 0, Z 2 = D. 

The asymptotic vdW interaction between large ob- 
jects is commonly found by summing the zero-point en- 
ergies of the weakly-damped plasmons of the combined 
system, minus that of the isolated systems 0, @. For 
T = OK, however, the isolated graphene layers do not ex- 
hibit weakly damped plasmons at small surface-parallel 
wavevector q. Therefore we need to go back to the more 
general and basic RPA correlation energy formulation, 
to which the sum-of-plasmon-energies approach is an ap- 
proximation, valid only when plasmon poles near the real 
frequency axis dominate (see e.g. ||J). The adiabatic 
connection-FDT gives the "cross-correlation" energy of 
the present well-separated system as (see e.g. 0): 



^ Jo ^ 



dX(J c (q,X,D)-J c (q,X,oo)) ,(9) 



/>OC />OC 

J c (q,X,D)=e 2 / du dzdz' exp(— q|z — z'\) x 

Jo J-oo 

\ X x{q, G = G' = 0, z, z', iu)} — [A — > 0]) (10) 



(1 



where x\ is the combined interacting response obtained 
above for two graphene planes separated by distance D. 
In the cross correlation energy 10 for distant graphenes, 
only the G = G' = component of the response func- 
tions are required because the coulomb interaction be- 
tween G 7^ fluctuations on the two layers is of order 



exp( 



G 



D) 



0(exp(-7rZ?/a)). This is << 1 when 
the layer separation D greatly exceeds the intraplanar 
spacing a between carbon atoms. ((TO and Q) give 



J c c r > RPA (q,X) 



-2V n (q) 



J c ,rpa(q,\D) - J c ,RPA(q\\,X, 00) 



2tt jo 



du(2V 12 (q)Xi2\{q\\,iu) 
(Xiia(A«>w) -Xiia(oo,3,m/)) (11) 



where Vu, V12 are defined in JSJl. Because every term 
in an expansion of (|llfl in powers of the small quantity 
V12 is at least of order = (9(exp(— 2qD)), we de- 
duce that the cross energy for distant slabs D — > 00 
is dominated by small parallel wavenumbers q of order 
D^ 1 << a -1 . For D >> awe can therefore approximate 
V11 , V12 and xo assuming that q is small compared with 
any G . In this limit (JSJ) is dominated by the G " = 
terms, Vu w 2ne 2 /q, V 12 « Vu exp(-qD). From 0, 



and (nj we tncn nave Jc r (q, A) = 2-£-F x {Dq), where by 
substituting u = v q sinhx and defining the dimension- 
less constant b = we 2 /Hvq = 12.14 5 for graphite we find 



F x (y) = - / 

Jo 



XbE+E_ 



1 



(c + XbE + ) (c + XbE_) 



Xb 



dx. 



Here c = cosh a; and E± = 1 ± exp(— y). Then the in- 
teraction energy per unit area is (in Gaussian esu units) 



E7 = / dXET = -Be 2 D~ 3 



(12) 



where B - (2tt) f dX J °° y 2 F x {y)dy = 2.003 6 x 
10 -2 is dimensionless and independent of D. 

Eq H12fl is a principal result of the present work. It 
shows that the asymptotic van der Waals attraction en- 
ergy, at T — OK, between parallel graphene planes dis- 
tant D, falls off like D~ 3 . This is to be compared with 
£)-5/2 f or 2D metallic planes [flU^, and Z?~ 4 for planes 
of finite-gap 2D insulators: the D~ 4 law follows equally 
from a sum of atom-atom R~ e contributions or from the 
zero-point energy of gapped plasmons. Clearly the gap- 
less graphene planes behave in this respect more like met- 
als than insulators, despite the lack of undamped 2D plas- 
mon modes on a graphene sheet. These examples illus- 
trate the non-universality of the vdW functional as well 
as its nonlocalitv|30j. Adding an exchange-correlation 
kernel f xc to the calculation of the response function in 
eq. © would modify the coefficient in Eq (|12fl but not 
the D~ 3 power law dependence. 

Another interesting case, possibly more accessible ex- 
perimentally, is the interaction between a graphene layer 
and a weakly metallic 2D layer. This could represent, 
for example, the outer part of the energy curve for inter- 
action between a layer being peeled off a solid graphite 
surface, and the last layer of the remaining solid, in an 
exfoliation experiment. Because of a weak overlap of or- 
bitals in the z direction perpendicular to the graphene 
planes, each layer of the solid has small pockets of elec- 
trons and holes and can reasonably be modelled by a 2D 
electron gas in the present context. 

We therefore applied the above RPA method to 
the non-overlapping interaction between an isolated 
graphene layer and a metallic graphite layer with Fermi 
energy e F = 0(0.02 eV). For D » Dq — , T(J , 
0(1 nm), the energy per unit area is (c.f. Q12[0 



2 v 2 _ 



e:: 



-Ce 2 D- 3 ln(D/D ) 



(13) 



Here C is a dimensionless constant. Like the case of two 
non-metallic graphene planes, the result l|13H disagrees 
with standard theories. 

In summary, our new calculations (see (|12(l . (|13|l . 30] ) 
reinforce our main point that the asymptotic (large-D) 
dispersion energy is a highly non-universal function of 



4 



separation D between non-overlapping subsystems. In 
particular we have highlighted the inadequacy, in this 
distant regime at least, of the usual sum of contribu- 
tions of form CqR~ 6 between pairs of atoms distant R, 
and therefore of various functionals that reduce to this 
asymptotic limit as R — » oo. A finite sum of multi- 
pole, or triplet and higher terms will also not reproduce 
what we have discussed. This inadequacy seems to occur 
when the component systems are (i) metallic (or have a 
zero electronic Bloch bandgap), (ii) of large spatial ex- 
tent in at least one dimension, so that long-wavelength 
(low-q) charge fluctuations are possible, and (iii) low- 
dimensional (of nanoscopic dimensions in another spatial 
direction), so that the electron-electron screening is re- 
duced compared to 3D bulk metallic systems, leaving a 
divergent polarizability at low frequency and wavcnum- 
ber. This physics is expected to be important in pre- 
diction of the cohesive energy of graphite, intercalated 
graphite, graphitic hydrogen storage systems, bundles 
of metallic nanotubes or nanowires, biomolecules and a 
large variety of weakly bound compounds ("soft mat- 
ter" ) . Although it is only in the widely-separated regime 
that the relevant long-wavelength electronic response is 
dominant, these low wavenumbers are still likely to be 
relevant near the overlapped equilibrium configuration. 
We have discussed some partial evidence for this conclu- 
sion, from existing jcllium layer calculations [12]. These 
used RPA-like nonlocal seamless microscopic theory that 
treats covalent, metallic and van der Waals effects on 
an equal footing. A similar level of theory is required 
for more realistic non-jellium models of soft matter, and 
this is leading to challenging numerical calculations [32| . 
Simplified versions [f| will probably have to take explicit 
account of large-scale geometry and/or nonlocal entities 
such as electronic bandgap. We note finally that our 
considerations might affect the analysis of some seminal 
experiments^^, E(l | concerning graphitic cohesion, be- 
cause these relied at some point on theory involving a 
sum of R~ e contributions. 
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